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Abstract. We present a semiclassical two-fluid model for an interacting Bose gas 
confined in an anisotropic harmonic trap and solve it in the experimentally relevant 
region for a spin-polarized gas of 87 Rb atoms, obtaining the temperature dependence 
of the internal energy and of the condensate fraction. Our results are in agreement 
with recent experimental observations by Ensher et al . 
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Bose-Einstein condensation (BEC) has recently been realized in dilute vapours of 
spin-polarized alkali atoms, using advanced techniques for cooling and trappingjl], ||, 
[|, ^ . These condensates consist of several thousands to several million atoms confined 
in a well which is generated from nonuniform magnetic fields. The confining potential 
is accurately harmonic along the three Cartesian directions and has cylindrical 
symmetry in most experimental setups. 

The determination of thermodynamic properties such as the condensate fraction 
and the internal energy as functions of temperature is at present of primary interest 
in the study of these condensates j|, p). The nature of BEC is fundamentally 
affected by the presence of the confining potential J6| and finite size corrections are 
appreciable, leading for instance to a reduction in the critical temperature 0, [j| |s|, |iofl. 
Interaction effects are very small in the normal phase but become significant with the 
condensation-induced density increase. The correction to the transition temperature 
due to interactions has been recently computed by Giorgini et al pj . 

The temperature dependence of the condensate fraction was recently measured || 
for a sample of around 40 87 Rb atoms, the observed lowering in transition 
temperature being in agreement with theoretical predictions within experimental 
resolution. In the same work the internal energy was measured during ballistic 
expansion and found to be significantly higher in the BEC phase than predicted by 
the ideal -gas model. While the increase is easily understood as a consequence of the 
interatomic repulsions, a quantitative estimate is still lacking. 

In this work we present a two-fluid mean-field model which is able to explain the 
above-mentioned effects, giving results in agreement with experiment for both the 
condensate fraction and the internal energy as functions of temperature. 
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We describe the condensate by means of the Gross-Pitaevskii (GP) equation for 
its wave function ^(r), 

?i 2 V 2 

*(r) + F ext (r)#(r) + 2gn 1 (rW(r) + g^ 3 (r) = uV(r) (1) 

2m 

where g = 4Trh 2 a/m, a being the scattering length, V QXt (r) — muj 2 (x 2 + y 2 + A 2 z 2 ) / 2 is 
the confining potential and n\{r) is the average non-condensed particle distribution. 
The factor 2 in the third term arises from exchange jl2| and we neglect the term 
involving the off-diagonal density of non-condensed particles. Following Bagnato et 
al Jl3| we treat the non-condensed particles as non-interacting bosons in an effective 
potential V (r) = V CVLt (r) + 2gni(r) + 2g^ 2 (r) . Thermal averages are computed with 
a standard semiclassical Bose-Einstein distribution in chemical equilibrium with the 
condensate, i.e. at the same chemical potential /i. In particular, the density ni(r) is 



m(r) 



1 f d 3 p 



( 2 ^) 3 J exp { (|i + V*(r) - /i) /k B T] - I 
= E - P {-^(r) - MksT} . (2) 

We fix the chemical potential from the total number of particles N 

where Nq = J ^ 2 {r)d 3 r and the semiclassical density of states is 



p(E) = ^^l JE-V<*(r)d 3 r. ,4) 



V° !t (r)<E 



This completes the self-consistent closure of the model. 

Equation (Q) can be solved analytically in the experimentally relevant situation 
N ~ I0 4 ^I0 5 and a/a ± ~ HT 2 , where a± = y/K/ rruo. Except for a small region close 
to the phase transition the interaction parameter Nga/a± entering the GP equation 
is large and the kinetic energy can be neglected. This yields 



>2 , . fi-V"* t (r)-2gn 1 (r) 



^2 {r) = w ^w 0(Ai _ F ext (r) _ 2 5ni ( r )) ( 5 ) 
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where 9{x) = (I) for x < (x > 0). The present strong-coupling solution neglects 
the condensate zero-point energy Eq = huo(\ + A/2). As Giorgini, Pitaevskii and 
Stringari|ll| pointed out finite-size effects are thereby excluded. 

Before presenting the complete numerical solution of the self-consistent model 
defined by equations (2)-(|J) and comparing its predictions with existing experimental 
data|^|, we display perturbative solutions at zero- and first-order. 

An approximate semi-analytical solution can be obtained by treating perturba- 
tively interactions involving the "dilute gas" of non-condensed particles. To zero order 
in gni(r) we have 
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and equation (Q) gives 
1 



Po(E) = 



TrX(hiu) 
for fx > and 



2^JI(E-pf 2 + £ 2 arctm/^ + [2^-E) 2 \n - ^~ 2/i ' 



(7) 



for p < 0. The self-consistent zero-order solution is then completed by equation (||). 
We remark that no assumption of weak interactions within the condensate has been 
made. 

We now proceed to compute the first order correction to the above zero-order 
solution. We take 

*»(r) = ^(0-2^(0) _ yext _ 

and expand equation (^) to first order in g[n\ (r) — n\ (0)] . The choice of the expansion 
parameter g[n\(r) — ni(0)] ensures that the perturbative expansion is regular, since 
the correction to p vanishes where the zero-order term vanishes. With the additional 
approximation g^ 2 {r) ~ pd(p) in the first-order term we get 



P {E) = p (E-2g ni (0)) + 5 Pl (E-2g ni (0) - p9(p)) (10) 



where 



k R T\ 3/2 ^ £; e J>e(-M)/fc fl r r / 3 ^ s 



Fi -.2, 



2' ' J fc B T 



(11) 



i-Fi being the Kummcr confluent hypcrgcometric function. The self-consistent first- 
order solution is then completed by equation (|^). 

We have solved numerically the simplified two-fluid model, first treating the 
parameter gni(0) to all orders (equations 2-||), then to zero order (equations [| and 
H ||) and finally to first order (equations || and P[ll]). Each case involves solving the 
integral equation @ to obtain p as a function of N and T; the non-perturbative 
solution also involves the local nonlinear problem posed by equations (2) and (g). 
The small differences between our three results justify a posteriori a perturbative 
treatment. Experimental parameters are taken from the work of Ensher et al Q: 
A = \/~8, N = 40000 and a/a± = 0.0062. We have verified numerically that our 
results depend weakly on N in the region explored in the experiments and therefore 
have used a fixed N = 40000 in all our computations. We use as energy units the 
semiclassical ideal-gas critical temperature fc^To = hui(NX/((3)) 1 ^ 3 , £(3) ~ 1.202 
being the Riemann zeta function. 

Figure [t] compares the temperature dependence of the condensate fraction Nq /N 
with the experimental results of Ensher et al g. Lowering of the transition 
temperature due to interactions is clearly visible, even if the smoothness of our results 
around the transition prevents a precise assessment of an interaction-induced shift in 
T c from the numerical solution. It should be noticed that the strong-coupling solution 
of the GP equation is not valid for T close to T c , since it requires No 3> a±/a ~ 160, 
and that our mean field model does not include critical fluctuations. Both effects 
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Figure 1. Condensate fraction obtained in the two— fluid model compared with 
the experimental data of Ensher et al [p| (diamonds) and with the ideal gas result 
(dotted curve) . We present results obtained from the zero-order solution (full curve) , 
from the first-order perturbative treatment (dashed curve) and from the complete 
numerical solution (circles). The inset is an enlargement of the region around T c . 



being relevant only in a narrow window around T c |Tl|, we expect our results to be 
meaningful in most of the temperature range. Recently Giorgini et al [ p| solved 
numerically the Popov approximation to the finite-temperature generalization of 
the GP equation within a semiclassical WKB approximation. Their results for the 
temperature dependence of the condensate fraction are in very good agreement with 
the predictions of our more naive model except for \T — T c \/T c < 0.05, where they find 
a sharp change in the slope of Nq(T). Their result for the interaction-induced shift 
in critical temperature ST C /T C ~ — \.33N 1 / 6 a/a± — —0.048 is also in good agreement 
with our curves. 

Figure || reports our results for the temperature dependence of the internal energy. 
We remark that the experimentally measured quantity is the sum of the kinetic energy 
and of the interaction energy, not including the confinement potential energy due 
to the rapid switching off of the trapping potential Jl4[. The average single particle 
energy (E) nc = J Ep(E)dE/ {exp [(E — /i)/igT] — 1} obtained from the semiclassical 
density of states contains twice the interaction energy, and - assuming that on average 
the kinetic and potential terms are equal - is twice the measured quantity. The 
kinetic energy of condensed atoms is negligible in our strong-coupling limit and 
their interaction energy per particle is (E) c = j9 J ^(rjcPr. The quantity directly 
comparable to the experimental data is therefore E — ((E) (N — Nq)/2 + (E) )/N, 
which we plot in figure || obtaining good agreement with the measured values. The 
calculated internal energy does not contain any sharp feature at transition, paralleling 
the result discussed above for the condensate fraction. Correspondingly the rapid rise 
in the specific heat is considerably smoothed with respect to the ideal-gas result (see 
figure ||). Apart from this small region around transition, our results on E(T) above 
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Figure 2. Sum of kinetic and interaction energy, as defined in the text, obtained 
in the two— fluid model compared with the experimental data of Ensher et al B 
(diamonds) and with the ideal gas result (dotted curve). We present results obtained 
from the zero-order solution (full curve), from the first-order perturbative treatment 
(dashed curve) and from the complete numerical solution (circles). The straight line 
is the classical Maxwell-Boltzmann result. The inset is an enlargement of the region 
around T c . 



and below T c imply a significant reduction of the increase in specific heat across the 
phase transition. 

In conclusion, we have presented a mean-field, semiclassical two-fluid model 
and discussed its perturbative and non-perturbative solution in the experimentally 
relevant parameter range. Our results on the temperature dependence of the 
condensate fraction and of the internal energy are in agreement with recent 
experimental measurements, accounting for the pronounced increase in internal energy 
with respect to the noninteracting boson case measured below T c . We have also verified 
that our model reproduces the results obtained for the condensate fraction with a more 
refined theory by Giorgini et al . 
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Figure 3. Specific heat C = dE/dT obtained in the two-fluid model compared with 
the ideal gas result (dotted curve). We present results obtained from the zero-order 
solution (full curve), from the first-order perturbative treatment (dashed curve) and 
from the complete numerical solution (circles). The straight line is the classical 
Maxwell-Boltzmann result. The anomalous behaviour obtained from the first-order 
perturbative solution near T c is an artifact. 
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